Nonthermal fixed points and solitons in a one-dimensional Bose gas 
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Single-particle momentum spectra for a dynamically evolving one-dimensional Bose gas are anal- 
ysed in the semi-classical wave limit. Representing one of the simplest correlation functions these 
give information about possible universal scaling behaviour. Motivated by the previously discovered 
connection between (quasi-)topological field configurations, strong wave turbulence, and nonthermal 
fixed points of quantum field dynamics, soliton formation is studied with respect to the appearance 
of transient power-law spectra. A random-soliton model is developed to describe the spectra ana- 
lytically, and the analogies and difference between the appearing power laws and those found in a 
field theory approach to strong wave turbulence are discussed. The results open a view on solitary 
wave dynamics from the point of view of critical phenomena far from thermal equilibrium and on a 
possibility to study this dynamics in experiment without the necessity of detecting solitons in situ. 

PACS numbers: 03.75.Lm 47.27.E-, 67.85.De 



<Z3 

c3 



i 

-a 
c 
o 
o 



in 

en 

en 

o 

CN 

> 

X 



I. INTRODUCTION 



Many-body systems far away from thermal equilibrium 
can show a much wider range of characteristics than equi- 
librium systems. However, they buy this abundance for 
the price of transiency. Among the wealth of possible 
non-equilibrium many-body configurations most interest- 
ing candidates for theoretical and experimental study are 
those at which generic time-evolutions get stuck for an 
extraordinarily long time. This is in particular possible 
at critical points where universal properties like power- 
law behaviour of correlation functions appear, leading 
to slowing-down phenomena as infrared modes become 
dominant. Out of equilibrium, fluid turbulence is among 
the earliest described as a critical phenomenon [1-4]. 

In the realm of quantum physics, far-from-equilibrium 
many-body dynamics, from the formation of Bosc- 
Einstein condensates in ultracold gases to quark-gluon 
plasmas produced in heavy-ion collisions and reheating 
after early-universe inflation, exhibits many interesting 
phenomena. In this context, increasing attention [5-18] 
is being given to wave turbulence phenomena [19-21]. 
Critical points far from equilibrium, so-called nonthermal 
fixed points, were proposed [6] and related to strong wave 
turbulence [7-10] and the formation of quasi-topological 
defects [15-17]. Such defects play an important role in 
superfluid turbulence, also referred to as quantum tur- 
bulence (QT), which has been the subject of extensive 
studies in the context of helium [22, 23] and dilute Bose 
gases [24-27]. In contrast to eddies in classical fluids, 
vorticity in a superfluid is quantised [28, 29], and the 
creation and annihilation processes of quantised vortices 
are distinctly different [22, 23]. 



'Electronic address: t.gasenzer@uni-heidelberg.de 



Considering a single vortex or vortex line, geometry 
imposes a power-law shape of the angle-averaged flow 
velocity away from the vortex core. This spatial power- 
law dependence implies the power-law momentum spec- 
trum predicted to occur at the nonthermal fixed point 
[15, 17]. In this article we focus on solitary waves in one- 
dimensional quantum gases and show that these exhibit 
scaling behaviour in analogy to the universal properties 
of vortex ensembles in two- and three-dimensional sys- 
tems. Self similarity here is due to the absence of a scale 
in a soliton which has zero width and is randomly posi- 
tioned in space. See Refs. [30, 31] for related studies in 
relativistic field theory. 

Solitons are non-dispersive wave solutions which can 
arise in many non-linear systems spanning a wide range 
from the earth's atmosphere [32] or water surface waves 
[33] to optics [34]. The characteristics of single or few 
solitons in ultracold Bose gases have been studied during 
the last decade regarding their movement and interaction 
in traps [35-38] , their formation and creation [39-42] and 
their decay [43, 44], see Ref. [45] for a review. Here, we 
propose to detect and characterise soliton excitations in 
ultracold gases by measuring the single-particle momen- 
tum spectra for such systems. We emphasise that these 
can give a strong indication for the presence of solitons 
even in cases where they cannot be observed in situ. 

Superfluid turbulence plays an important role in the 
context of the kinetics of condensation and the devel- 
opment of long-range order in a dilute Bose gas. This, 
as well as turbulence in its acoustic excitations has 
been discussed in Refs. [46-50], and more recently in 
Refs. [15, 17, 18, 51]. A possible observation of QT in 
ultracold atomic gases presently poses an exciting task 
for experiments [52-54]. We stress that the experimen- 
tal study of superfluid turbulence and, more generally, of 
nonthermal fixed points in ultracold Bose gases is central 
for the understanding of the processes important during 
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FIG. 1: Snapshots of a single run of the non-linear classi- 
cal field equation, showing solitons which oscillate inside a 
trapped one-dimensional ultracold Bose gas. The gas is ini- 
tially noninteracting and thermalised, with T = 36Ck<Jho, in 
a trap with oscillator length Iho = 8.5 (in grid units). At 
time t = the interaction is switched to gm — 7.3 x 10~ 3 , 
and a cooling period using a high-energy knife is applied, 
see Sect. Ill for details on grid units, the chosen param- 
eters and protocol. The panels show the one-dimensional 
colour-encoded density distribution as a function of time. 
Top left panel: The initially imposed interaction quench 
causes strong breathing-like oscillations and the creation of 
many solitons. Top right panel: Breathing oscillations have 
damped out, leaving a dipolar oscillation of the bulk distri- 
bution in the trap. Clearly distinct solitons have formed. 
Bottom left panel: After the end of the cooling period (here 
at t = t c — 9.1 x 10 3 ) only a few solitons are left oscillat- 
ing within the oscillating bulk. Bottom right panel: A single 
soliton is left at late times. 



the build-up of coherence and degeneracy in particular 
also in one-dimensional systems [55-57]. 

In the following we study the formation of soliton exci- 
tations in trapped one-dimensional Bose gases by means 
of simulations in the classical- wave limit of the underlying 
quantum field theory. We analyse the momentum spec- 
tra of the time-evolving system and discuss them with 
respect to nonthermal fixed points discussed in field the- 
ory. In Section II, we develop a model of independent, 
randomly positioned grey solitons, locally being solutions 
of the Gross-Pitaevskii equation describing a degenerate, 
one-dimensional dilute Bose gas. Single-particle momen- 
tum spectra are derived for such systems, both in a homo- 
geneous system and within a trapping potential. A pro- 
tocol for the formation of such soliton ensembles within a 
trapped gas is described in Section III, and the resulting 
states are analysed with respect to the predicted momen- 
tum spectra and power-law signatures of a nonthermal 
fixed points. Our conclusions are drawn in Section IV. 



II. MOMENTUM SPECTRA OF SOLITON 

ENSEMBLES 

In the classical- wave limit, a dilute ultracold Bose gas 
is well described by a positive definite Wigner phase- 



space distribution function W[<j>, 4>*] for the complex field 
4> at each point in position or momentum space [58, 59]. 
The dynamics of the gas is determined by the time evo- 
lution of this Wigner function according to the classical 
field equation 
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which is identical in form to the Gross-Pitaevskii equa- 
tion (GPE) for the field expectation value. Here, m is the 
mass of the bosons, V an external trapping potential, and 
(/id the coupling constant in one spatial dimension (ID). 
The ID Gross-Pitaevskii equation (1) possesses quasi- 
topological soliton solutions which may travel with a 
fixed velocity but are non-dispersive, i.e., stationary in 
shape [45]. For positive coupling constant g 1D > soli- 
tons are characterised by an exponentially localised den- 
sity depression with respect to the surrounding bulk mat- 
ter and a corresponding shift in the phase angle ip of the 
complex field <j> = |^>| exp{iip}. Depending on the depth 
of this depression, the soliton is either called grey or, for 
maximum depression, black. On the background of a 
homogeneous bulk density n it is described by 



4> v {x,t) = y/n 



1 



tanh 
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vt is the position of the soliton at 

-1/2 



where X s (t) — x 

time t. Here, £ = [2mngiu]^ 1/2 is the healing length, 
and we express time in units of the inverse speed of sound 
c s = [ngi-a/m] 1 ' 2 , i.e., t = i/c s , and drop the overbars. 



is the 'Lorentz factor' corresponding to 



7 = l/vT 

the velocity v of the grey soliton in units of the speed of 
sound, v = v/c a = \4> v (vt, t)\/y/n. Being related to the 
density minimum, v is also termed the 'greyness' of the 
soliton, ranging between (black soliton, \4> v (vt,t)\ — 
v^/n = 0) and 1 (no soliton, \<f> u (vt,t)\ — \/n). 

In Sect. Ill we will study the formation and evolution 
of soliton excitations in trapped one-dimensional Bose 
gases by means of the GPE (1). An example of one run 
is shown in Fig. 1: A sudden initial quench of the cou- 
pling 51D creates strong oscillations of the bulk gas in the 
trap and lets solitons form out of the short-wave-length 
collective oscillations. These solitons are seen as black 
lines surviving within the oscillating gas for a long time. 

In this article we are mainly interested in the char- 
acterisation of the ensemble of solitons emerging in our 
simulations in terms of single-particle spectra in momen- 
tum space. Hence we refer to Sect. Ill for more details on 
the simulation protocol and results and start here with a 
discussion of the possible spectra. 

The single-particle momentum spectrum is determined 
by evaluating ensemble averages 



f»(M) = MM)| S 



'ensemble 



(3) 



over a large number of runs. Fig. 2 shows n(k, t) at three 
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FIG. 2: Momentum spectrum n(k, t) at different times t 
as indicated, with statistical errors (ensemble averages over 
1000 runs). Momentum scales defined by the inverse heal- 
ing length, fc^, as well as by the density of solitons, fc„ s , are 
marked by dashed lines. A power-law dependence ~ k~ 2 is 
indicated by the black solid line. Here, a larger trap oscilla- 
tor length was chosen, lh = 17, such that around the outer 
turning point of the breathing oscillation the solitons become 
diluted and thus the scales k na and k^ sufficiently far away 
from each other to allow for a k~ 2 power law to be seen in 
between (red points). Momentum k and time t are measured 
is in grid units as defined in Sect. Ill where also all other 
simulation parameters are given. 



different times during the period where the initial breath- 
ing oscillations are still present. Note the double loga- 
rithmic scale. At low momenta, the spectrum shows a 
plateau while at high momenta it falls of exponentially. 
At the point of time corresponding to the red curve, the 
cloud is further expanded such that solitons are more 
separated from each other which causes an intermediate 
power law dependence n(k) ~ fc~ 2 to appear, as indi- 
cated by the straight line. The plateau, the power-law 
and the exponential decay form the characteristic signa- 
ture for solitons which we discuss in more detail in the 
following. 



A. Random-soliton model: uniform gas 

To obtain an analytical understanding of the possible 
spectra in the context of nonthermal fixed points we dis- 
cuss, in the following, the case of a dilute ensemble of 
well-separated solitons with random velocities and po- 
sitions. The wave function of a single grey soliton in 
a homogeneous bulk background condensate is given in 
Eq. (2). Using this we write down an expression for a 
set of N s uncorrelated solitons with density minima at 
{xi,i = 1, ...,N S } and (dimensionless) velocities z/;, on a 
background of constant bulk density n: 



^\x,t) = V^U [n-^^ix-x,) 



(4) 



Note that due to the neglection of correlations this field 
does in general not represent a solution of the GPE in 
which the solitons remain non-dispersive. 

We make use of the assumption that the ensemble is 
dilute, i.e., that the distance between each pair of neigh- 
bouring solitons is much larger than the healing length. 
This assumption, for grey solitons, is not valid as soon 
as two oppositely moving solitons encounter each other 
but for simplicity we will assume that these collisional 
configurations can be neglected in view of a majority of 
well-separated solitons. Since, for any i, 4>' Vi (x — Xi) = 
d$ Vi (x — Xi)/dx ~ as \x — Xi\ 3> 1 we can rewrite the 
spatial derivative of the field (4) as 



J& 



E 



5(x - »«(*)) JJe*^*-*'^ 



0,0). (5) 



Here Xi(t) = Xi(0) — Vit, /3j = 2 arccos za,- , * denotes the 
convolution over the spatial dependence on x, 0{x) the 
Heaviside function, and we have neglected an irrelevant 
overall phase. Note that the sign of (3i indicates the di- 
rection of the propagation of the ith soliton. The term 
in square brackets in the second line of Eq. (5) is pro- 
portional to the spatial derivative of the field describing 
an ensemble of N s infinitely thin solitons (£ —> 0) , at the 
positions {x^, 



(iv s 



'(M)-E 






(6) 



where n s is the number of solitons per unit length and 
the prefactor containing 7, takes into account that the 
phase jump by exjp{if3 i 6(x — Xi)} is itself proportional 
to a theta function Ji6(x — Xi). Note that although the 

derivative <J>c_£q (Xi,t) gives a sum of terms, each being 
proportional to a delta distribution, only one of these 
remains when evaluated at Xi, which gives the term in 
square brackets in Eq. (5). We take the Fourier transform 

of (0^^ (x) 0( JVs ) (y)) with respect to x — y, integrate 
over R= (x + y)/2, and divide by k 2 , 
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Here, P is the probability for finding a soliton with grey- 
ness v = cos(/3/2), and averaging over the random po- 
sitions of all solitons other than those at Xi and Xj has 
been done in order to obtain the exponential decay of the 
coherence function. Combining Eq. (7) with the Fourier 



transform of <j> v (x, 0) 

4> v (k,t) = i 



27m 



2iru6(k) 



V2< 



sinh (n^fk^/y/2) 



(8) 



one derives the single-particle momentum distribution for 
a set of N s solitons defined by greyness and position, 
{vi,Xi \i = 1,, 



n(k,t) 



n 
4n? 



N s }, as 
dk 1 
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4k' n s Rea 
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(9) 



Here, the inverse volume 6(0) = L^ 1 appears as we first 
choose the (/>*. and <p v . in Eq. (7) different and take the 
identity limit only at the end. a is the average over all 
an = (1 - exp{0i})/2. 

Assuming the dependence of %/ sinh(7r7jfc£/v2) on Vi 
to be negligible we obtain an approximate stationary dis- 
tribution 



i(k) 
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, (10) 



with a yet to be determined parameter 7. For black soli- 
tons (vi = 0) one obtains the exact expression 



mu = 



4n s n (irk£,) 2 /2 
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fc2 sinh 2 (irk£/V2 



For an ensemble of grey solitons of identical |i^| = ^, trav- 
eling with probabilities P into the positive x-direction 
and Q = 1 — P into the negative direction one finds 



n(k) 



An*n 



4n 2 7 - 4 + [k + 2(1 - 2Q)n s ^7- 1 ] 2 
(nk0 2 /2 



sinh' 



(wjkt/y/2 



(12) 



Here, 7 = (1 — ^ 2 ) -1 / 2 . To demonstrate the applicability 
of the above analytic expressions we construct ensembles 
of phase-space distributions of spatially well-separated 
solitons in a box with periodic boundary conditions and 
compute the ensemble average (3). These simulations are 
done on a ID grid of TV = 2 14 sites, generating 5 x 10 3 
configurations for taking ensemble averages. For this we 
multiply N s = 20 single-soliton solutions (2) with posi- 
tions Xi and greyness i>i chosen randomly according to a 
given phase-space distribution. To make sure that their 
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FIG. 3: Single-particle momentum spectrum as a function 
of fc on a double-logarithmic scale for an ensemble of 5 x 
lO 3 configurations with N s — 20 solitons each distributed 
according to a flat distribution across phase-space defined by 
the positions in the box and the maximum greyness |^ m ax| = 
0.99. We chose £ = 8a a . Solid (black) squares: numerical 
ensemble averages, solid (black) line: Eq. (10) with a = 0.7, 
7 = 1.05, shown up to fc£ ~ 1.3 where a finite-size effect sets 
in due to the random set of solitons not matching periodic 
boundary conditions. For comparison, results for iV s = 20 
purely black solitons (red squares and line) are shown as well 
as for N a — 20 solitons with fixed greyness \u\ = 0.707, i.e., 
7 = 1.4, (blue squares and line) choosing an equal number 
of right- and left-movers, P — Q = 1/2. The comparison 
validates the approximate expressions (10)-(12) which exhibit 
scaling behaviour in the regime of momenta fc ns 7~ 2 "C k <§C 
fc 57 -\ 



relative distance on the average is much larger than their 
widths we chose the phase-space distribution to allow for 
a maximum greyness, \i/i\ < l^maxl < E such that the 
diluteness criterium requires an approximate minimum 



(11) box length of L = 4(1 



1-1/2AT 



Fig. 3a shows the single-particle momentum spectrum 
n(k) on a double-logarithmic scale for an ensemble of 
5 x 10 3 configurations with N s — 20 solitons each dis- 
tributed according to a flat distribution across phase- 
space defined by the positions in the box and the maxi- 
mum greyness j^maxl = 0.99. Solid (black) squares repre- 
sent the results of the numerical ensemble average while 
the solid line corresponds to the analytical formula (10), 
with fitted parameters a = 0.7, 7 = 1.05. Compare 
this to the analytical average a = 2/3. For compari- 
son, we give the results for the same number of purely 
black solitons (red squares and line) as well as for a fixed 
greyness \u\ = 0.707 (blue squares and line), choosing 
an equal number of right- and left-movers. The com- 
parison validates the approximate expressions (10)-(12) 
which exhibit scaling behaviour in the regime of momenta 
kn s l~ 2 < A:< fc C 7" 1 . 

In Fig. 4 we show the single-particle momentum spec- 
trum for an ensemble of 5 x 10 3 configurations with N s — 
20 solitons each distributed according to a phase-space 
distribution with an unequal weight for solitons with 
positive v (right-movers) and negative v (left-movers). 
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FIG. 4: Single-particle momentum spectrum as a function of 
k on a double-logarithmic scale for an ensemble of N = 5 x 10 3 
configurations with N 3 = 20 solitons each distributed ac- 
cording to a phase-space distribution with an unequal weight 
for solitons with positive v (right-movers) and negative v 
(left-movers). The greyness is uniformly distributed within 
= fmin < v < t'max = 0.99. All other parameters are chosen 
as for Fig. 3. Red squares and black line show the ensemble 
average and compare with Eq. (12) for Q = 0.2, 7 = 1.179, 
and v — 0.412, shown up to fc£ ~ 1 where a finite-size ef- 
fect sets in due to the random set of solitons not matching 
periodic boundary conditions. For comparison, results for 
N s — 20 solitons with fixed greyness v = 0.707, i.e., 7 = 1.41, 
(blue squares and line) are shown. 



We specifically restricted the greyness to the interval 
= z^min < v < f max = 0.99 and besides that chose 
the same parameters as for Fig. 3. Black squares and line 
show the ensemble average and compare with Eq. (12) for 
Q = 0.2, 7 = 1.179, and v = 0.412. For comparison, re- 
sults for N„ = 20 solitons with fixed greyness v — 0.707, 
i.e., 7 = 1.41, (blue squares and line) are shown. 



B. Relation to nonthermal fixed points and 
vortical excitations in superfluids 

Universal power-law behaviour in a many-body system 
far from equilibrium points to the appearance of turbu- 
lence phenomena. In the following, we summarise a few 
basic ideas of wave-turbulence theory for quantum gases 
and discuss the power-law spectra derived for soliton en- 
sembles in this context. 

A dilute, degenerate Bose gas is compressible such that 
collective sound-wave excitations can occur. This allows 
for so-called weak wave-turbulence which can occur in 
a regime where kinetic theory applies: The wave-kinetic 
equation for a Bose gas, which describes the evolution 
of the system under the collisional interactions between 
different wave modes, is known to have non-trivial sta- 
tionary solutions. These solutions are nonthermal and 
exhibit a power-law dependence of mode occupations on 
momentum [19-21]. As in fluid turbulence, such solu- 
tions imply that energy flows from, e.g., small to large 



momentum scales. In between these momentum scales of 
the source and the sink, the energy passes through the 
so-called inertial interval, where the distribution over mo- 
menta is stationary and follows a power law with a uni- 
versal exponent predicted by weak-wave-turbulence the- 
ory [19-21]. The stationary solution is metastable, i.e., 
requires a constant and equal flow in and out of the in- 
ertial interval. We note that weak wave turbulence does 
not bear, in general, vortical excitations. It is defined by 
scaling and stationary transport between different scales 
and hence can also occur in one spatial dimension where 
vortices are absent. 

Applying the above to a degenerate Bose gas one faces, 
however, the problem that the description in terms of 
wave-kinetic equations breaks down in the infrared (IR) 
regime of long wavelengths where amplitudes, i.e., talk- 
ing about a Bose gas, single-particle occupation numbers 
grow large and the description in terms of, e.g., elastic 
two-to-two collisions becomes unreliable [19]. As a con- 
sequence of this, so-called strong wave turbulence is ex- 
pected to occur in the IR regime. Recent developments 
presented in Refs. [6-10, 13] allow one to set up a uni- 
fying description of scaling, both in the ultraviolet (UV) 
regime, where the wave-kinetic and Quantum Boltzmann 
equations apply, and in the IR limit. Recall that at a 
critical point the IR modes dominate the system's be- 
haviour. In this IR regime, new scaling laws were found 
by analysing non-perturbative Kadanoff-Baym dynamic 
equations, as in the UV, with respect to nonthermal sta- 
tionary power-law solutions [6-8, 10]. These solutions 
were termed nonthermal fixed points. Analogous predic- 
tions for dilute Bose gases were given in [9], proposing 
what was termed strong matter-wave turbulence in the 
regime of long-range excitations. 

In [15, 17], it was then shown that these nonthermal 
fixed points can also be understood, in two and three 
dimensions, in terms of vortex excitations of the super- 
fluid: In the infrared limit of large wave numbers the 
incompressible, superfluid component of the gas domi- 
nates, and the predicted IR power laws appear due to 
the algebraic radial decay of the flow- velocity around the 
vortex cores. As a consequence, within a window in mo- 
mentum space, which is limited by the inverse mean dis- 
tance between different vortices and the inverse core size, 
the single-particle occupation number spectrum shows 
the power law predicted in [9] . 

Before we discuss this in more detail let us first turn 
back to the soliton spectra derived in the previous sec- 
tion. Assuming an equal number of solitons traveling 
with positive and negative velocities, P = Q = 1/2, 
i.e., assuming Imo = 0, the single-particle spectrum 
(12) is characterised by a maximum of two scales. Con- 
sider the case n s <C 7 / '{\/2-k^) . For momenta greater 
than the reduced soliton density but smaller than the 
reduced inverse healing length, fc„ s 7 -2 <C k <C fc^7 _1 , 
with fc„ s = 2n s , and k^ — \/2/(iv£), the momentum dis- 
tribution exhibits a power-law behaviour, n(k) ~ k~ 2 . 
This reflects, first, the random position of the kink-like 



phase jump across the center of the soliton, and second, 
that these momenta cannot resolve the spatial width of 
the kink. In other words, looking within a spatial win- 
dow of size between j 2 k7 and 7&" 1 , the appearance of 
a single sharp solitonic phase jump inside the window is 
observed in a random manner. Hence, for any of these 
window sizes, the system looks identical, it appears self- 
similar. This self-similarity is at the base of the scaling 
momentum distribution. Already for a single soliton, the 
momentum distribution does not know anything about 
the position of the kink (this information appears as a 
phase in the momentum-space Bose field which is irrele- 
vant for n(k)), and thus the single-soliton distribution is 
self-similar, too. 

For black solitons, the self-similar scaling region is lim- 
ited by the scales k ng and k^. Below k na , the distribution 
is constant because too low wave numbers cannot resolve 
the kink-structure. This corresponds to the first-order 
coherence function decaying exponentially in space, with 
the decay scale set by the soliton distance l/n s , 



'(x)cj>(y)) ~ exp{-2n s |a; - y\}. 



(13) 



Above fcg, the momentum spectrum resolves the finite 
width of the soliton density dip which results in an ex- 
ponential suppression of the mode occupations. We re- 
call that also in equilibrium, at sufficiently high temper- 
atures where quasiparticle mode occupations are large, 
n< qp )(fc) = k B T/(c s k) > 1/2, the first-order coher- 
ence function decays exponentially, g^'(s) ~ exp[— s/ro], 
where the scale is set by the coherence length r = 
2n/(mfcsT). Hence, the corresponding momentum spec- 
trum has the same shape n cq (k) ~ 2r ~ 1 /(r ~ 2 + k 2 ) as for 
a set of random thin solitons. We emphasise, however, 
that the transition scale fc„ s above which scaling ~ fc~ 2 
sets in (see, e.g. (11)) can be made larger than for a ther- 
mal ensemble by increasing the soliton density n s above 
the inverse thermal coherence length 1/tq- This allows 
to identify non-equilibrium soliton vs. thermal scaling in 
experiment. 

We finally compare the universal and non-universal as- 
pects of the soliton momentum spectra found here with 
the corresponding spectra in d = 2 and 3 dimensions. As 
discussed in detail in Refs. [15, 17], the universal expo- 
nent C = d + 2, d = 2,3, found for the particle spectra 
n(k) ~ fc~^ during the dynamical relaxation of an ini- 
tially strongly quenched gas reflected the appearance of 
vortices in two and vortex lines in three dimensions. In 
two dimensions, this can be seen in an easy way looking 
at the flow velocity field v ~ V</? = eg/r at the distance 
r from the core of a singly quantised vortex, where eg 
is the local tangential unit vector and tp the phase angle 
of the complex Bose field. The r-dependence implies a 
/c _1 scaling of |v(fc)| and thus a k~ 2 scaling of the kinetic 
energy E(k) ~ k 2 n(k) ~ v(k) 2 , i.e., n(k) ~ fc~ 4 [17, 60]. 
Similar arguments lead to n(k) ~ k~ 5 for the radial mo- 
mentum distribution in the presence of a vortex line in 
three dimensions [17]. Extending these arguments, the 



scaling n(k) ~ k~ d ~ 2 was shown to appear for ensembles 
of randomly positioned vortices/vortex lines in a range of 
momenta l~ x < k < £ _1 between the inverse of the inter- 
vortex distance l v and the inverse of the healing length £ 
which is a measure for the core width. 

As pointed out above, the power-law spectrum n(k) ~ 
k~ , in turn, had been predicted by use of nonpertur- 
bative field-theory methods in [6, 9] where it resulted for 
a strong-wave-turbulence cascade in the IR, characteris- 
ing the scaling behaviour at a nonthermal fixed point. 
This cascade was shown in [17] to be caused by particles 
being transported towards the IR where they build up 
high mode occupations and thus coherence in the sample, 
see also [46-51]. Note that in the picture of the evolving 
Bose field this momentum-space transport corresponds 
to the mutual annihilation of vortices and anti-vortices 
in the system which results in an increase of the inter- 
vortex distance and thus of the range over which phase 
coherence is established. 

Having recalled all this, we note that there is a discrep- 
ancy between the predicted scaling ~ k~ d ~ 2 , which was 
found consistent with the vortex picture for d — 2 and 3, 
and the scaling ~ k~ 2 obtained here for the solitons in 
d = 1. To gain more insight into this issue we consider the 
spatial decay of the phase coherence for a system in two 
dimensions over distances considerably larger than the 
intervortex spacing. In analogy to the soliton ensembles 
discussed in Sect. II C, we consider randomly positioned 
and well-separated vortices. One finds that the decay of 
the phase coherence follows the same exponential law, 



(x)0(y)) ~ exp{-2n v ,i |x - y|}, 



(14) 



where n v ,i is the one-dimensional uniform vortex density 
along the straight line through x and y. Fourier trans- 
forming (14) with respect to x — y results in a momentum 
spectrum n c (k) ~ (4n 2 x + fc 2 ) -3 / 2 which scales as fc~ 3 
for momenta considerably larger than the inverse vortex 
distance n v ,i = l/l v - Analogously, one finds n c (k) ~ fc~ 4 
for vortex line tangles in three dimensions. We use the 
subscript 'c' to distinguish these spectra from the n(k) 
discussed above. 

The apparent contradiction between the respective 
scalings of n(k) and n c (k) whose exponents differ by 1, is 
resolved by observing the following: While the exponen- 
tial decay (14) is valid over large distances |x — y| 3> l v , 
it is the algebraic dependence of the flow field v as a 
function of distance from the vortex core which matters 
on length scales below the inter-vortex distance, caus- 
ing the steeper power law n(k) ~ k~ d ~ 2 at momenta 
k ^$> 1/iv Hence, we find an important qualitative dif- 
ference between the physical properties underlying the 
universal scaling n(k) ~ fc~ 2 in one dimension and the 
scalings n(k) ~ fc~ d ~ 2 in d = 2 and 3 dimensions: While 
black solitons in one dimensions are at rest and no parti- 
cle flow can occur, in higher dimensions, transverse flow 
circling around the vortex cores gives rise to an addi- 
tional contribution to the kinetic energy. As far as this 



transverse flow dominates over a possible additional lon- 
gitudinal flow component for which v • Vv 7^ 0, see [17], 
this comparison also holds when allowing for grey soli- 
tons, which are moving opposite to the (longitudinal) 
particle flow across the soliton dip. 

In summary, there is a principal difference between the 
scalings of the momentum spectra in one- and in higher 
dimensions, giving rise to a deviation of the scaling ~ fc~ 2 
in d = 1 from the field-theory prediction k~ d ~ 2 which, 
in turn, is valid in d = 2 and 3 dimensions. To recover 
this discrepancy within a field-theory approach to strong 
wave turbulence is beyond the scope of this article. 
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C. Random-soliton model: trapped gas 

In our dynamical simulations we will consider soliton 
formation in a harmonic trapping potential rather than in 
a homogeneous system, see Sect. III. We therefore need 
to take into account, in the random-soliton model, the 
inhomogeneous bulk distribution of the gas. 

Assuming a sufficiently shallow harmonic potential, 
W£ = (muj^ 2 )- 1 / 2 = (2g 1D nuj ho )- 1 / 2 > nj 1 we 
can describe the Bose field in local-density approxima- 
tion with respect to a bulk density distribution given in 
Thomas-Fermi approximation, n^p(x) ~ n Q [l — (x/R) 2 ], 

with R — v2c s /(whoC) = 2giDn/who being the Thomas- 
Fermi radius in units of £. We take the maximum den- 
sity no large enough to ensure k^ ^> fc„ s for solitons not 
too close to the edge of the cloud. In such a bulk den- 
sity distribution, single solitons oscillate harmonically be- 
tween classical turning points where the solitons "touch 
ground", i.e., momentarily turn black [44]. Their oscilla- 
tion frequency is by a factor of ~ l/y2 smaller than the 
trap frequency, w s ~ u>\ io /y/2. 

In leading order in e ~ x s /x s , the field of a single 
soliton can locally be written in the simple form given in 
Eq. (2), with v, and thus 7 replaced by local quantities, 



v -» v(x s ) = v s 



7 -> j(x s ) = [1 - v(x s ) 2 



-1/2 



(15) 



evaluated at the posi- 
m ax is the maxi- 



tion x s = x s (t) of the soliton [44]. v s 
mum greyness the soliton acquired in the center of the 
trap, and a; S)max = ^ S ,max-R is the distance of the soliton's 
turning point from the trap center. Only solitons whose 
velocity does not exceed the Landau critical velocity, i.e., 
for which ^ S]max < 1, can oscillate in the trap for more 
than a quarter of the period T s = 2it/uj s . This limits the 
maximum greyness at a distance x from the trap center 
to a range between and Vm ax (x) = 1 — (x/R) 2 . 

At a given time t, we assume a particular set {xi(t}} of 
N s well-separated solitons across the trapped gas. The 
single-particle momentum spectrum corresponding to an 
ensemble of such sets depends on the distribution of the 
solitons over the greyness for each position in the trap. 



FIG. 5: Single-particle momentum spectrum n(k) on a 
double-logarithmic scale for an ensemble of 5 x 10 3 config- 
urations with N s — 20 solitons each distributed according to 
a flat distribution in phase-space {x/R, v}, circularly sym- 
metric around (x = 0, v = 0) with radius R s — 1. We chose 
£ = 4a s , iVo. Solid (red) squares: results of numerical en- 
semble average; solid line through these points: analytical 
spectrum (19), with 7 = 1.1. (Blue) squares and line: Cor- 
responding results for the same number of black solitons dis- 
tributed randomly across the trap. 



This distribution is best visualised in phase space which 
is parametrised by the (x,v), or, equivalently and in di- 
mensionless form, by (x/R, v), with both, x/R and v 
ranging between —1 and 1. In this space, the trajectory 
of a single soliton is a circle with radius v s . max which is 
traced out with constant angular velocity u> s . Hence, a 
stationary distribution of the N s solitons is given by a 
circularly symmetric distribution in phase space, i.e., a 
distribution over the different possible maximum grey- 
nesses z/ Sjmax or turning points x Sjmax . The simplest as- 
sumption would be that of a uniform distribution of the 
solitons in phase space, which amounts to a uniform dis- 
tribution over the different possible v at each distance 
x from the trap center and an integrated soliton density 
distribution n s (x) ex 1 — (x/R) 2 . 



Following the above considerations we can obtain ap- 
proximate expressions for the momentum spectrum. For 
instance, for a uniform density P[x, v\ = n s ^/R s of soli- 
tons within a radius R s = R s /R in phase space, i.e., 
for all (x,v) with y/x 2 + v 2 < i? s < 1, with n Sj o = 
N s /(R s tt), the first-order coherence function for thin soli- 
tons becomes 

<*<"•> (a) V"-%)>«-o 

= \/ n TF(x)n TF (y) exp{-2n Si0 / dza(z)} (16) 



with a local average dephasing of 



<\ (,)•) = _ _ / dv \\ — v ±wyl-y 2 



2fl, 



-•/S^-x 2 




(17) 



Using this, the integral over a to linear order in a: reads 



f*dza(z) = a«x + 0(x 3 ) 



,(i) 



1 - i?2/3. This ap- 



proximation is best for R s — > 1, in which limit we obtain 

~ v/n TF (x)n TF (y)e- 4 ^l s -«l/( 3 ^). (18) 

Analogously, one calculates the exponential dephasing 
factor for more complicated soliton phase-space distribu- 
tions. Taking the above results together one derives the 
momentum distribution in local-density-approximation 
as the convolution of the spectrum for a homogeneous 
distribution of thin solitons with the Fourier transform 
of the bulk density, multiplied with the momentum spec- 
trum of a single soliton: 



n(k) ~ btf(^) * 



4n sfl a 



(i) 



(7T 7 fc0 2 /2 



4<o c 



,(1)2 



sinh 



TT-yk£/y/2) 
(19) 



where the parameter 7 is to be determined and 7k- denotes 
the convolution with respect to k. 

In the case that single soliton distributions contribut- 
ing to the ensemble are not uniform throughout phase 
space, the ensemble-averaged n(k) would rather be a sum 
of (k — > — fc)-asymmetric distributions such that on the 
average the momentum distribution can have local max- 
ima at finite \k\ > 0. 

To study the quality of the above analytic expressions 
we construct ensembles of phase-space distributions of 
spatially well-separated solitons inside a harmonic trap 
and compute the ensemble average (3) . For this we mul- 
tiply N s single-soliton solutions (2) with positions Xi and 
greyness i>i chosen according to a given phase-space prob- 
ability distribution and ensure that their relative distance 
on the average is much larger than their widths. 

Fig. 5 shows the single-particle momentum spectrum 
n(k) on a double-logarithmic scale for an ensemble of 
5 x 10 3 configurations with N s = 20 solitons each dis- 
tributed according to a flat distribution in phase-space 
{x/R, v}, circularly symmetric around (x = 0, v = 0) 
with radius R s = 1. Solid (red) squares represent the 
results of the numerical ensemble average while the solid 
line corresponds to the analytical formula (19), with 
7 = 1.1. For comparison, we give corresponding results 
for the same number of black solitons distributed ran- 
domly across the trap (blue squares and line). 



III. SOLITON SPECTRA IN DYNAMICAL 
SIMULATIONS 

A. Soliton formation and tracking in position space 

In the following we study the formation of soliton en- 
sembles by use of semiclassical simulations, with Gaus- 
sian noise for the initial field modes. Moments of the 
phase-space probability distribution at a later time are 
determined by sampling the initial distribution, propa- 
gating each realisation according to the GPE, and av- 
eraging over many such trajectories [58, 59]. At the 
initial time, we take the gas to be noninteracting and 
thermalised and impose an interaction quench. To al- 
low the emerging collective excitations to form solitons 
at a desired density we furthermore apply evaporative 
cooling by opening the trapping potential at the edges 
in a controlled fashion. During the first, cooling pe- 
riod, t < t c , the potential is given by the inverted 
Gaussian V(x,t) = mujl U{t){l-exp[-x 2 /2U(t)}} with 
its maximum being ramped down by sweeping U(t) — 
Uq + (U c — Uo)t/t c linearly in time from Uq to U c . At 
the same time, highly energetic particles near the edge 
of the potential are removed by adding a loss term 
iT(x,t)/2 = iT 00 [V(x,t)/U(t)] r /2 to the trapping po- 
tential. Thereafter, during the interval t c < t < i max 
the loss is switched off and the potential is ramped up 
again to harmonic shape accross the extension of the 
gas, U(t) = U c + (IU - U c )(t - t c )/(t max ~ t c ). We 
choose r = 10, U = 2.75, T^ = 0.1 U , U c = C/ /3, and 
U max = 10 Uq. The times t c and t max vary and are given 
in the following. This protocol corresponds to the one 
used in Ref. [61]. Different cooling schemes have been 
used in experiments, see, e.g., [55, 57], but as we are pri- 
marily interested in the one-dimensional dynamics, we 
here restrict ourselves to purely ID calculations. 

For the simulations, we map the system onto a grid 
of N — 1024 lattice sites with a lattice constant a s . If 
not stated otherwise, quantities are given in grid units 
based on o s , and the parameters chosen are a dimension- 
less coupling constant g 1D = 2ma s gm — 7.3 x 10 -3 , a 
cooling time t c = t c /(2ma 2 ) = 9.1 x 10 3 , and a harmonic 
oscillator length Zho = ^ho/a = 8.5. Lattice momenta are 
k n = 2sin(Trn/N), n G {-N/2,..., N/2 - 1}. We drop 
overbars in the following. 

Three stages of the induced dynamical evolution can 
be observed, see also Fig. 1: 

Initial Oscillations. See the top left panel of Fig. 1. 
Following the interaction quench potential energy is 
transferred to kinetic energy. One observes strong 
breathing-like oscillations of the gas. These oscillations 
decay on a timescale of t rj 2 x 10 3 , leaving a dipolar 
oscillation of the bulk in the harmonic trap. Solitons are 
formed in the wake of the decaying breathing oscillations. 

Solitonic Regime. See the top right panel of Fig. 1. The 
initial collective oscillations have largely decayed except 
for an overall dipole mode, and many solitons appear. 
The bottom left graph shows the evolution around t = 
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FIG. 6: Snapshots of a single run of the non-linear classi- 
cal field equation, showing solitons which oscillate inside the 
trap, thereby showing signs of mutual scattering and passing 
through each other (examples marked by letters). Parameters 
are chosen as in Fig. 1. 



t c = 9.1 x 10 3 when only very few solitons have survived. 
The solitons oscillate in the trap, being nearly black at 
the edges and grey in the center of the trap corresponding 
to a nonzero velocity. Fig. 6 magnifies a short period of 
the evolution. On mutual encounters, the solitons get 
phase-shifted, such that collisions show signs of scattering 
or passing through each other. Collisions with different 
such shifts are marked by letters A and B in Fig. 6. 

Final stage: At times t ^> t c , a soliton is still visible, 
see the bottom right panel of Fig. 1. Comparing runs we 
find different numbers of solitons remaining during the 
late stage. 

The smallest time scale is the oscillation period in the 
trap Xho ~ 230, which leads to an initial collective os- 
cillation with period T osc in at ion ~ 300 (cf. Fig. 1). The 
collective breathing motion dies out after r sa 2000. The 
oscillation period of a soliton in the Thomas Fermi bulk is 
T s sa \/2T\ m fa 320. The longest time scale in our setup is 
the cooling time t c = 9100. Comparing these time scales 
to the total time of the simulation, at the end of which 
solitons are still present, we see that the solitons are 
quasi-stationary in the system. They emerge soon after 
the initial quench and remain throughout the whole evo- 
lution while thermalisation of the high-momentum modes 
is proceeding as we will see in the following. 



B. Number of solitons after cooling ends 

In order to study the statistics of the solitons emerging 
during the evolution we have set up an efficient tracking 
algorithm which identifies the trajectories of the solitons 
oscillating in the gas. The algorithm scans the wave 
function for density minima coinciding with a phase 
jump around them. Fig. 7 shows the evolution of the 
mean number of solitons, for an ensemble of 200 runs. 
The three stages described above can be identified. The 
strong initial oscillations give an oscillating number of 
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FIG. 7: Time evolution of the mean number N s of solitons, 
for an ensemble of 100 runs. Strong initial oscillations occur 
during the breathing like bulk oscillations after the interaction 
quench, until t ~ 10 3 . The number of solitons decreases while 
the gas is evaporatively cooled. After the end of the cooling, 
at t — t c the decay is considerably slowed down and the num- 
ber of solitons remains largely stable. Three different cooling 
times and two ramp speeds are shown, with t c — 3.9 x 10 3 
(red), t c = 4.55 x 10 3 (blue), and t c = 9.1 x 10 3 (black). 



solitons until t ~ 10 3 . The number of solitons decreases 
while the gas is evaporatively cooled. After the end of 
the cooling, at t = t c the decay is considerably slowed 
down and the number of solitons remains largely stable. 
Three different cooling times and two ramp speeds are 
shown, with t c = 3.9 x 10 3 (red), t c = 4.55 x 10 3 (blue), 
and t c = 9.1 x 10 3 (black), where the same speed is 
chosen to obtain the blue and black data. 



Kibble and Zurek have predicted that the number of 
defects created in the near-adiabatic crossing of a phase 
transition scales with the crossing rate according to a 
power law which depends on the universal properties of 
the transition [62, 63]. This was studied numerically in 
[61] using the cooling protocol described above. While 
the interacting gas was chosen to be in thermal equilib- 
rium initially, with a temperature well above the critical 
point, we start our simulations, motivated by earlier work 
on vortex dynamics [15, 17], with an interaction quench 
driving the system strongly out of equilibrium. To com- 
pare the dynamics induced in this way with the results 
of [61] we show, in Fig. 8, the dependence of the number 
of solitons created on the cooling ramp time t c . We find 
that, within the error bars which indicate the variance 
over 200 runs, the data is rather fitted by an exponential 
dependence N s (t c ) — /q exp{— jt c } than by a power law 
N s (t c ) = go/t c as predicted in Ref. [39]. We emphasise 
however that in our system, solitons mainly form during 
the initial stage following the interaction quench. 
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FIG. 8: The number of solitons N s (t c ) at the end of the cool- 
ing period, t = t c , as a function of t c , for ensembles of 200 runs 
each (black dots and error bars). The function /o exp{ — jt c }, 
with /o = 9.7 and 7 = 9.77 • 10~ 5 , was fitted with \ 2 = 0-006 
(black line). The function got^ 1 , with go — 1.55 • 10 5 , was 
fitted with x 2 — 1-28 (red dashed line). Hence, other than for 
the Kibble-Zurek scheme of Ref. [61], cooling after the initial 
quench results in an exponential dependence of N 3 (t c ) on t c . 



C. Time evolution of single-particle spectra 

We finally discuss the relaxation dynamics with re- 
spect to the evolution of the respective single-particle 
momentum spectra (3). The initial state chosen in the 
simulations is given by a thermal canonical ensemble of 
distributions over the single-particle eigenstates of the 
trap. In Fig. 9 we show the momentum spectrum at time 
t = 5- 10 3 . Solitons have formed at high density such that 
the scales k Us ~ 0.15 and k^ ~ 0.2 are close together, as 
indicated in the graph. The solid line represents a fit of 
the analytical model spectrum (12). A Rayleigh- Jeans 
tail is absent as the cooling is still on. Due to the prox- 
imity of k nB and k^ no fc~ 2 power law is seen in between 
the low-energy plateau and the high-energy exponential 
fall-off. 

In Sect. II, Fig. 2 we showed the spectrum for a wider 
trap with l\ 1Q — 17 which allows the solitons to be di- 
luted more across the trap and results in fc„ s ~ 0.15 and 
&£ ~ 0.65. This allows for a k~ 2 scaling to appear clearly, 
indicating a self-similar random distribution of the soli- 
tons. For t > t c , the gas enters the final stage: The num- 
ber of particles and the energy are now conserved and the 
healing length is 'frozen out'. Fig. 10 shows the develop- 
ment of n(k, t) from t > t c to late times. Once the cooling 
and thus the removal of particles with high energy is ter- 
minated, a transport process from low to high momenta 
starts and thermalisation takes place. The influence of 
solitons is still present with two effects: First, there are 
still solitons in the gas for the times displayed in Fig. 10 
which contribute with their spectral profile to the total 
spectrum and broaden the plateau at low momenta up 
to fcj ~ 0.1. Second, the momentum distribution starts 
to oscillate between situations with a stronger weight on 
the positive and on the negative side. Fig. 11 shows how 
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FIG. 9: Momentum spectrum n(k, t) at time t — 5 • 10 3 with 
statistical errors (ensemble average over 100 runs), k and t 
are in grid units as defined in Sect. Ill A where also all other 
simulation parameters are given. Solitons have formed at high 
density such that the scales fc„ s and k^ are close together, 
as indicated in the graph. Solid line: fit of the analytical 
model spectrum (19), with n s ,o = 0.076, 7 = 1, a~l = 0.036 
where Gaussian of width a was used to describe the bulk 
distribution in position space. A Rayleigh-Jeans tail is absent 
as the cooling is still on. Due to the proximity of k„ s and k^ 
no k~ 2 power law is seen in between the low-energy plateau 
and the high-energy exponential fall-off. 
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FIG. 10: The spectrum n(fc) is shown for three different times 
from t « t c = 9 x 10 3 to t = 10 6 , for an ensemble of 5 x 10 3 
(10 3 for t = 10 6 ) realisations, k and t are in grid units as de- 
fined in Sect. Ill A where also all other simulation parameters 
are given. Momentum scales defined by the inverse Thomas- 
Fermi radius R as well as by the healing length £ are marked 
by dashed lines. A power-law dependence ~ k~ 2 of the ther- 
mal tail is indicated by the black solid line. 



the remaining solitons, which oscillate in the trap, influ- 
ence the spectrum with their own momentum appearing 
as the oscillating maxima in the spectrum. 
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FIG. 11: Time evolution of the spectrum shown on a linear 
color scale at late time 4.95 < 10 -5 £ < 5.0, for a single run. 
The oscillatory pattern appears despite a large number of runs 
contributing to the statistical ensemble. A small number of 
runs with a few oscillating solitons dominate the averaged 
spectrum. 



IV. CONCLUSIONS 

We have studied the formation of dark solitary waves 
in one-dimensional Bose-Einstein condensates as well as 
their relaxation dynamics towards equilibrium. The cor- 
responding single-particle momentum spectra were pre- 
dicted in the framework of an instantaneous model of 
well-separated grey solitons, whose width is considerably 
smaller than their mutual distances in the bulk. For com- 
parison with these predictions, semiclassical simulations 
of the relaxation dynamics of one-dimensional Bose gases 
after an initial interaction quench and a cooling period 
were used to determine the respective spectra numeri- 
cally. The so found spectra compared well with the an- 
alytical predictions, giving insight into the many-body 
dynamics from the point of view of universal properties 
and critical physics far from equilibrium. We emphasise 
that the particular protocol used to produce the solitons 
is irrelevant in that the properties characterising the fixed 
point do not depend on how it is reached. 

We have discussed the power-law behaviour of the 
momentum spectra, which appears in a range of mo- 
menta between the inverse of the inter-soliton distance 
and the inverse healing length, with regard to the uni- 



versal scaling laws predicted in non-perturbative field- 
theory approaches to strong wave turbulence. In the one- 
dimensional case studied here, the derived power-law ex- 
ponent C = 2 differs by 1 from the exponent C, — d+2 pre- 
dicted for a system in d spatial dimensions and previously 
recovered in the scaling due to vortex excitations in a 
d = 2 and 3-dimensional Bose gas. We trace this discrep- 
ancy back to the different flow patterns possible in d > 1 
versus d = 1 dimensions: While in the one-dimensional 
gas, particle flow cannot choose its orientation except for 
a sign, vortical excitations in two and three dimensions 
are characterised by flow circling around the vortex cores. 
This flow is transverse, i.e., it changes its strength per- 
pendicular to its direction. It dominates, at sufficiently 
low energies and momenta, over any additional longitudi- 
nal flow caused by compressible sound-wave excitations. 
Moreover, for geometric reasons, the transverse-flow ve- 
locity field changes algebraically in space and causes the 
single-particle momentum spectra to scale as predicted 
in strong-wave-turbulence theory. 

We point out that the single-particle momentum spec- 
tra discussed here could be used in experiment to study 
solitary-wave dynamics in one-dimensional Bose gases 
without the necessity to detect solitons in situ. Study- 
ing in this way universal properties during the relaxation 
dynamics from a non-equilibrium initial state or under a 
constant driving force opens a new access to strong wave 
turbulence and nonthermal fixed points. 
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